/**
* @Arquivo: lagrangian.c
* 
*	Geração e Solução do problema de coloração usando a formulção de representantes 
*
* @author: Marcio Costa Santos
* @data: 20/05/2010
*
*/

#include "lagrangian.h"
#include "pl.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

extern int n;//o numero de vertices
extern int m;//o numero de arestas
extern int k;//o k da coloração
extern long va[];//a ordem dos verties

CPXENVptr env[MAX_SIZE];//o ambiente do problema
CPXLPptr prob[MAX_SIZE];//os problemas de cada vertice
double x[MAX_SIZE][MAX_SIZE];//o vetor das soluções de cada problema
double bestLB;//o melhor limite inferior encontrado
double bestUB;//o melhor limite superior encontrado
double lb;//o limite inferior corrente
double ub;//o limite superior corrente
double pi;//passo do lagrangeano
int noI;//melhoria
int handle;//auxiliar
double coef[MAX_SIZE];//os coeficientes da função objetivo
long map[MAX_SIZE*MAX_SIZE]; //mapeia as variaveis nos seus mu

/*Função de ordenação para o qsort*/
int compar(const void* a, const void* b){
	int v = *(int*)a;
	int u = *(int*)b;
	return mu[v] - mu[u];
}

int define(){
	int i,status;
	pi = 2.0;
	noI = 0;
	lb = ub =  0;
	bestLB = (-1)*n;
	bestUB = n;
	for(i = 0; i < n; i++){
		mu[i] = Lx[i] = 2.0;
		env[i] = CPXopenCPLEX(&status);
	}
	if(status!= 0){
		printf("problemas para abrir cplex");
		fflush(NULL);
	}
	return initPL();
}


int makePL(long v){
	return addConstrainsForVertex(v, &(env[v]) ,&(prob[v]));
}

double calculateUB(){
	int V[n];
	int cor[n];
	int i,j;
	int color = 1;
	int adj[n];
	int cand = -1;
	for(i = 0; i < n; i++){
		V[i] = i;
		cor[i] = 0;
		adj[n] = 0;
	}
	qsort(V, n, sizeof(int), compar);
	cor[0] = 0;
	for(i = 1; i< n; i++){
		for(j = 0; j < i; j++){
			if(!hasEdge(V[i], V[j])){
				adj[cor[j]] = 1;
			}
		}
		for(j = color - 1; j >= 0; j--){
			if(adj[j] == 0)cand = j;
		}
		if(cand == -1){
			cand = color;
			color++;
		}
		cor[i] = cand;
		cand = -1;
		for(j = 0; j < n; j++){
			adj[j] = 0;
		}
	}
	if(bestUB > color){
		bestUB = color;
	}
	return color;
}

int noImproviment(){
	noI++;
	if(noI > 30){
		pi = pi/2.0;
		noI = 0;
	}
	return 0;
}


double calculateLB(){
	int i;
	double value = 0.0;
	for(i = 0; i < n; i++){
		value = value + Lx[va[i]] + mu[va[i]];
	}
	if(value > bestLB){
		printf("valor do limite melhorado para %f \n", value);
		bestLB = value;
	}
	else noImproviment();
	return value;
}

	
int updateMuLambda(){
	double g[MAX_SIZE];
	double theta = 0.0;
	double theta_div = 0.0;
	int i,j;
	for(i = 0; i < n; i++){
		g[va[i]] = 1;
		for(j = 0; j <= i; j++){
			g[va[i]] -= x[va[j]][va[i]];
		}
		theta_div += g[va[i]]*g[va[i]];
	}
	if(theta_div >=  0.001)theta = pi*(ub - lb)/theta_div;
	else return 1;
	//calcula subgradientes
	for(i = 0; i < n; i++){
		mu[va[i]] = mu[va[i]] + theta*g[va[i]];
		if(mu[va[i]] < 0.0000){
			mu[va[i]] = 0.0;
			if(g[va[i]] <= 0.0) g[va[i]] = 0.0;
		}
	}
	theta_div = 0.0;
	//recalcula subgradente
	/*for(i = 0; i < n; i++){
		theta_div += g[va[i]]*g[va[i]];
	}
	if(theta_div >=  0.001)theta = pi*((1.05)*ub - lb)/theta_div;
	else theta = pi/n;
	//recalcula subgradientes
	for(i = 0; i < n; i++){
		mu[va[i]] += theta*g[va[i]];
		if(mu[va[i]] < 0.0000){
			mu[va[i]] = 0.0;
		}
	}*/
	return 0;
}

double solvePL(int v){
	int i = 0;
	Lx[va[v]] = 1 - solve(env[v], prob[v], v) - mu[va[v]];
	getX(env[v], prob[v], x[va[v]], v);
	if(Lx[va[v]] > 0 ){
		Lx[va[v]] = 0;
		for(i = 0; i < n; i++){
			x[va[v]][i] = 0.0;
		}
	}
	return Lx[va[v]];
}

int addObjFunc(long v){
	int i, j;
	j=0;
	for(i = v+1; i < n; i++){
		if(hasEdge(va[i], va[v])){
			coef[j++] = mu[va[i]];
		}
	}
	return changeObj(env[v],prob[v], coef, v);
}


double lagrangian(double* sol){
	int i;
	for(i = 0; i < n; i++){
		makePL(i);
	}
	while(pi > 0.005){
		for(i = 0; i < n; i++){
			addObjFunc(i);
			solvePL(i);
		}
		lb = calculateLB();
		ub = calculateUB();
		if(updateMuLambda() != 0)return bestLB ;
	}
	*sol = bestUB;
	return bestLB;
}


